SANEPIC: A Map-Making Method for Timestream Data From 

Large Arrays 



G. Patanchon, 1 ' 2 ^ P. A. R. Ade, 3 J. J. Bock, 4 ' 5 E. L. Chapin, 1 M. J. Devlin, 6 S. Dicker, 6 
M. Griffin, 3 J. 0. Gundersen, 7 M. Halpern, 1 P. C. Hargrave, 3 D. H. Hughes, 8 J. Klein, 6 
G. Marsden, 1 P. G. Martin, 9 ' 10 P. Mauskopf, 3 C. B. Netterfield, 10 ' 11 L. Olmi, 12 ' 13 
E. Pascale, 11 M. Rex, 6 D. Scott, 1 C. Semisch, 6 M. D. P. Truch, 14 C. Tucker, 3 
G. S. Tucker, 14 M. P. Viero, 10 D. V. Wiebe 11 



ABSTRACT 

We describe a map-making method which we have developed for the Balloon-borne Large 
Aperture Submillimetcr Telescope (BLAST) experiment, but which should have general applica- 
tion to data from other submillimeter arrays. Our method uses a Maximum Likelihood based ap- 
proach, with several approximations, which allows images to be constructed using large amounts 
of data with fairly modest computer memory and processing requirements. This new approach, 
Signal And Noise Estimation Procedure Including Correlations (SANEPIC), builds upon several 
previous methods, but focuses specifically on the regime where there is a large number of detec- 
tors sampling the same map of the sky, and explicitly allowing for the the possibility of strong 
correlations between the detector timestreams. We provide real and simulated examples of how 
well this method performs compared with more simplistic map-makers based on filtering. We 
discuss two separate implementations of SANEPIC: a brute- force approach, in which the inverse 
pixel-pixel covariance matrix is computed; and an iterative approach, which is much more effi- 
cient for large maps. SANEPIC has been successfully used to produce maps using data from the 
2005 BLAST flight. 

Subject headings: methods: data analysis — techniques: image processing — submillimeter — balloons 
1. Introduction 
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ogy differing dramatically between different sub- 
fields of astronomy. The method adopted depends 
on the form in which the data are gathered, and on 
the dominant source of systematic effects. From 
the optical to the near-IR one talks about combin- 
ing "frames" , along with measurements of "darks" 
and "flat-fields" . For the reduction of Cosmic Mi- 
crowave Background (CMB) data the now conven- 
tional method is to start from the principle that 
there is a linear algebra approach to solving the 
Maximum Likelihood problem. However, this has 
only been feasible up until now because of the 
limited number of detectors in the typical CMB 
experiment, and the fact that correlated signals 
among the detectors can be effectively ignored. 
Because of the rapid development of large bolome- 
ter arrays, the question which arises is: how does 
one adapt the CMB approach to dealing with sub- 
stantial numbers of detectors, and where there are 
significant cross-correlations of noise between the 
detector timestreams. 

Data from the Balloon-borne Large Aper- 
ture Submillimeter Telescope (BLAST, Devlin 
et al. 2004) represent a new challenge for bolo- 
metric timestream-to-map algorithms. Recent 
CMB experiments which use detectors similar to 
those used on BL AST, such as B OOMERANG 



I Cri ll et al.l 12003) and Archeops ( Benoit et al 
l2002h . only use a handful of separate bolometers. 
Furthermore, these experiments' off-axis designs 
lead to small correlations between detectors. Con- 
sequently, the correlations could be ignored at the 
map-making stage, and each detector timestream 
could be treated as an independent sub-set of the 
data. This has changed for BLAST, which has up 
to 139 detectors per band, with significant corre- 
lations induced by the on-axis design, as well as 
the higher frequencies of the observations. Just 
by itself, the large number of channels increases 
the impact of even small time stream correlations, 
the contribution from which does not integrate 
down with increasing number of detectors, unlike 
the uncorrelated noise. The high level of correla- 
tion (largely induced by temperature drifts in the 
obscuring secondary supports in BLAST) make it 
important that the correlations be handled care- 
fully in the map-making process. 

This paper describes a Signal And Noise Esti- 
mation Procedure Including Correlations (SANEPIC) 
which has been developed for the analysis of 



BLAST. This algorithm will also have applica- 
tion to many next generation experiments which 
will involve both noise correlations between chan- 
nels (including correlations from the atmosphere) 
and very large numbers of detectors. This includes 
the next generation of larger- format arrays for use 
in ground- and balloon-based instruments at mi- 
crowave and millimeter wavelengths, which will 
have typically thousands of detectors. 

This paper is arranged as follows. We next de- 
scribe the pertinent aspects of BLAST. In Sec- 
tion 3, the longest section, we set out our basic 
map-making method. Simulations we use for test- 
ing the method are described in Section 4, with 
results presented in Section 5, demonstrating the 
benefit of accounting for the correlated noise in 
BLAST- like observations. Finally some of the 
maps obtained from the June 2005 BLAST flight 
are presented in Section 6. 

2. BLAST observations of the submilli- 
meter sky 

The map-making procedure presented in this 
paper has been used to analyze the data from 
the BLAST 2005 flight. We will use BLAST as a 
specific example for the application of SANEPIC 
throughout the paper. We describe BLAST in Sec- 
tion 12.11 we then summarize the pre-processing of 
the data prior to map-making in Section 12.21 In 
Section [2~3l we derive a model of the data that we 
will use for the map-making. 

2.1. BLAST instrument and observations 

The Balloon-borne Large-Aperture Submilli- 
meter Telescope incorporates a 2-m primary mir- 
ror, and large format bolometer arrays operat- 
ing at 250, 350, and 500 /xm, designed to have 
144, 96 and 48 bolometers, respectively (of which 
139, 88 and 43, respectively, were used). The in- 
strument is described in detail in IPascale et al 



(|2007l ). The low atmospheric opacity at the oper- 
ating altitude of ~ 38 km allows BLAST to map 
the sky very quickly compared to ground-based 
experiments and to conduct large area shallow 
surveys as well as very deep surveys of the sky 
(|Devlin et al.ll2004) . The BLAST wavelengths are 
near the peak of the spectral energy distribution 
of cold galactic dust, which gives BLAST the abil- 
ity to conduct unique extragalactic and Galactic 
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submillimeter surveys with high spatial resolution 
and sensitivity. BLAST thus enables studies of 
the distribution of very high-redshift galaxies and 
of star forming regions in our Galaxy. 

The typical observing strategy consists of scan- 
ning the telescope back and forth in azimuth, cov- 
ering the entire field by slowly varying the ele- 
vation. Cross-linking of the data is assured by 
scanning the same field at another time of the 
day. Typica l scan ning strategies are given in 



Pascale et al.l (|2007l) . 

The first scientific flight of BLAST took place in 
June 2005 from the Esrange Arctic base in Swe- 
den to the Canadian Arctic. A total of ~ 100 
hours of data were taken in a variety of Galactic 
fields. They include a star fo rming region (Vulpec - 



ula) over 4 deg , described in lChapin et alJ (|20071 ). 



three other fields of similar size in the Galactic 
Plane (which will be the focus of future papers), an 
integration towards the ELAIS-N1 field (see Oliver 
et al. 2000), the Cas A supernovae remnant over 
about 0.5 deg 2 (Hargrave et al. in preparation), 
and sev eral compact Gala ctic and extra-galactic 
sources ( Truch et al.1 12007 ). Hereafter we refer to 
these as the BLAST05 data, to distinguish them 
from the data taken during the December 2006 
Antarctic flight. 

2.2. Time-ordered data pre-processing 

The processing of BLAST data from detector 
timestreams to the final map product involves sev- 
eral steps prior to map-making. Each of these 
steps is designed to remove a particular (or sev- 
eral) artifact (s) from the data, and sometimes re- 
quires iterating, since some effects need to be re- 
moved simultaneously. In the following, we sum- 
marize the main processing stages leading to the 
time-ordered segments which are used as inputs 
for the map-making process. 

We start by identifying events in the data which 
are sharply localized in time, such as spikes from 
cosmic rays hits and other spurious sources. We 
use a method which allows us to discriminate be- 
tween the different events depending on their sig- 
nature in the data. Spikes which involve only a 
single sample are flagged and the corrupted sam- 
ples are replaced by the average value of the sam- 
ples in the vicinity. The data are deconvolved 
from the low-pass filter applied by the readout 



electronics. This filter has a frequency cut off 
of approximately 35 Hz and is designed to avoid 
high frequency noise aliasing. The deconvolution 
is performed in Fourier space. In addition, we 
have applied a low-pass cosine filter which limits 
the noise power from increasing at very high fre- 
quency (above 38 Hz) due to deconvolution. We 
have checked that the noise power spectrum is rel- 
atively flat after these deconvolution operations. 
Finally, cosmic ray hits and other localized arti- 
facts in the data timestreams are detected and cut 
out. In order to avoid biasing our data products 
by having systematically more false event detec- 
tions located where the sky is bright (e.g., when 
scanning a point source), we iterate this process - 
we make maps starting from data which has been 
cleaned using the process described above, sub- 
tract the maps from each original data-set, and 
reprocess the data. The maps calculated at this 
intermediate stage are obtained by simply rebin- 
ning the data into pixels after strongly high-pass 
filtering. The filter applied is a Butterworth filter 
with a frequency cut off of 0.5 Hz, which is of the 
order of the knee frequency of the noise. Even if 
this operation suppresses most of the intermedi- 
ate to large scales of the sky signal in the maps, it 
does not very much alter the signal from localized 
sources, at least to the level of accuracy needed at 
this stage, and it has the advantage of removing 
most of the stripes in the maps due to 1// noise. 
We have verified that the resulting bias for bright 
calibrators is less than 1%. 

About 2% of the data from the BLAST05 flight 
was removed due to cosmic ray events. Most of 
the events affected a single detector timestream, 
although some events affected a whole array at 
the same time. Detected spikes (from cosmic rays 
but also other spurious effects) are flagged over 
small time intervals of typically 1 second in the 
data. One second is too large an interval to sim- 
ply ignore, so the corrupted data need to be re- 
placed with random noise generated in a way that 
as much as possible preserves the statistical prop- 
erties of the data. In the later map-making stage, 
which is partly performed in Fourier space, we as- 
sume continuity of the data. We could perform 
this gap-filling with a constrained noise realization 
(see e.g., Stompor et al. 2002). However, since 
the gap intervals are significantly smaller than the 
inverse of the knee frequency of the noise power 
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spectrum (0.3 Hz), the noise can be well approx- 
imated by the sum of a white component plus a 
straight line of some slope across the gap. Specif- 
ically, we generate white noise in each gap with 
a standard deviation measured from the data in 
the vicinity of the gap, and add a baseline with 
the parameters fit using 20 samples on each side. 
The white noise generation is for restoring as best 
as possible the stationarity of the data (generated 
samples are not reprojected to the map at the 
end). 

After having filled the gaps in the timcstrcams, 
we filter out very low frequency drifts which are 
poorly accounted for in the map-making proce- 
dure. A fifth-order polynomial is fit to the data 
and removed from each data segment in order to 
reduce fluctuations on timescales larger than the 
length of the considered segment which, depend- 
ing on the specific case, varies from 30 minutes to a 
few hours. These fluctuations are poorly described 
because of the limited number of Fourier modes, 
and would cause leakage at all timescales (for in- 
stance a gradient in the timestream is described by 
a wide range of Fourier modes), degrading the effi- 
ciency of map-making. Note that we experimented 
with various polynomials and other effective high- 
pass filters; we found that the results were not 
very sensitive to precisely how this is done, but 
some such filtering is certainly required. The de- 
gree of the polynomial was chosen empirically as a 
compromise between suppressing the artifacts and 
keeping the large scale signals in the final maps. 
Using simulations, we have checked that the ef- 
fect on the transfer function of the signal in the 
final map is weak. The resulting data segments 
are then corrected for the time- varying calibration 
(see Truch et al. 2007) u sing measurements of an 
internal calibration lamp (jPascale et al ]|2007i) . Fi- 
nally, the data segments are apodized at the edges 
over ~ 2, 000 samples and are high-pass filtered at 
5 x 10~ 3 Hz with a Butterworth filter. This filter- 
ing has very little effect on the final maps, since 
modes in the data at lower frequencies that are 
cut in this way are not expected to contribute sig- 
nificantly to the signals. We discuss the choices of 
filters further in Section l3~2l 

Accurate pointing reconstruction is a compli- 
cated procedure for balloon-borne telescopes, and 
this affects the map-making task through the 
pointing matrix (defined in Section 12. 3p . The 



pointing reconstruction proc edure is described in 
detail in lPascale et al. ( 2007 ). The next important 
step in reducing the BLAST data, which is cali- 



bration of the detectors, is detailed in Truch et al 



120071) 



2.3. Model of the data 

Having performed the cleaning procedure de- 
scribed in the previous sub-section, the resulting 
data timestreams can be modeled very accurately 
as the sum of pure signal and pure noise contribu- 
tions. The data for detector i observing at a given 
wavelength and at time sample t can be written 
as 

du = [Ai] tp s p + n it , (1) 

where p labels the pixels in the final map, Ai is the 
pointing matrix for bolometer i (whose elements, 
indexed with time t and pixel p, give the weight of 
the contribution of pixel p to the sample at time t 
for bolometer i), s p is the signal amplitude at pixel 
p, and the noise amplitude at time t for bolome- 
ter i is Ha. Summation over repeated indices is 
assumed here. 

Ideally, the clement [^4i]t p of the pointing ma- 
trix is equal to the value of the beam response 
b(R(r — ro)), where r points to the pixel p loca- 
tion, r~o is the location of the beam center at time 
i, and R is a rotation matrix which depends on 
the rotation angle at time t between the telescope 
and sky coordinate systems. In principle one could 
then recover a map of the sky deconvolved with the 
instrumental Point Spread Function (PSF). How- 
ever, in practice this would be unacceptably noisy, 
as well as computationally intractable, because of 
the prohibitive volume of data. Even although Ai 
might have mostly zero elements, it is nevertheless 
a huge matrix. It may be feasible to deconvolve a 
non-trivial beam response (e.g. like the BLAST05 
beams, as shown in Truch et al. 2007) at the same 
time performing the map-making step, through an 
approximate treatment of the sparse pointing ma- 
trix. But we did not pursue that approach here. 

In the simple case where the beam is symmet- 
ric, the map-making problem becomes tractable, 
provided one restricts oneself to reconstructing a 
map of the instrument-convolved sky. We can then 
consider s p in Equation [1] as the map of the sky 
convolved with the beam, and consequently Ai in- 
dicates simply where the detector points in the 
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sky at a given time. In this case, Ai is an ulti- 
mately sparse matrix with, in the BLAST case, 
simply a 1 in a single entry of each row. This 
approach, which has been conventional for CMB 
map-making (although with some adaptation for 
chopped data), is what we use in the following 
analysis. It gives no loss of information provided 
that the map pixels are sufficiently small, and one 
can simply assign all the flux from a bolometer's 
to the map pixel to which it points at each time 
interval. The requirement for accuracy is that the 
pixel size is smaller by a factor three or more than 
the Full Width at Half Maximum (FWHM) of the 
instrumental PSF; in that case the additional con- 
volution with the pixel shape gives negligible loss 
of angular resolution. 

The noise term nn in the model represents the 
sum of all contributions to the timestreams which 
do not reproject on the sky. This will in general 
include instrumental noise, fluctuations in atmo- 
spheric emission and other loading and unrecog- 
nized cosmic ray hits. In general, some of those 
noise contributions will induce strong correlations 
between detector timestreams. In this paper, we 
adopt a very general model of the noise where the 
noise covariance matrix: 



= (nu ■ n\, v ) 



(2) 



(for bolometer indices i, i' and time indices t, t') 
has possibly non-zero elements even for i i' . 
A key assumption, as we will see later, is that 
the noise is Gaussian (so that Nam' is sufficient 
to describe all the statistical properties of the 
noise) and stationary (constraining Na'tt 1 , see Sec- 
tions GTO and O . 

In the specific case of BLAST05 observations, a 
very significant correlation of the noise is found in 
the timestreams, and we have shown that a more 
constraining model provides a very good descrip- 
tion of the data. An independent component anal- 



ysis ( Delabrouille. Cardoso. &; Patanchonl2003f ) of 
the data enabled us to find that the noise and 
its correlations can be described to a high de- 
gree of accuracy by a noise component which is 
not correlated between detectors, together with a 
single common-mode component seen by all the 
detectors at a given wavelength (some correla- 
tions is also seen between detectors from differ- 
ent wavelength but we have chosen to treat each 
wavelength independently). The common part of 



the noise is instantaneous, meaning that the same 
common-mode noise is seen at the same time by 
all the detectors. In our model, the noise term in 
Equation [1] is then decomposed as: 



n lt = n it + o>iCt, 



(3) 



where the first term is the noise which is uncorre- 
cted between detectors and the second term rep- 
resents the common-mode component of the noise, 
rescaled by an amplitude parameter on which de- 
pends on the detector but not on time. This model 
can be generalized easily to deal with multiple 
noise components in timestreams. 

In the following section, we present a method to 
reconstruct s p given the data, in the framework of 
the linear model (Equation [T]) and in the presence 
of correlated noise between detector timestreams. 

3. Map-making method 

3.1. Maximum Likelihood map-making 

The use of Maximum Likelihood map-making 
techniques has been developed by many authors 
for applica tion to Cosmic Microwave Back g round 
data-sets (IWright et all Il996t iTeemarkl Il997t 
i Borrilll Il9~99l IPrunet et all IgOOd: iTegmark et al ' 



2000 : iFerreira fc Jaffd I2OOOI: iDore et all 12 001 



Natqli et alj200lt|Prunet et alj200lHDupac fc Giard 



20021: IStompor et al l 120021 : Iffinshaw et all 12003: 



Yvon & Mavetjl2005h . Some other approach 



arc 



more speci fic to destriping for Planck- like scanning 
strategies (iDelabrouillel 19981: iMaino et all [2002; 



2005 



2006; 



Keihanen et al.l 120051: Ide Gasperis et al.l 
Macias-Perez et al.l l200^ ~ Poutanen et al.l 
Ashdown et al .1 l2007h . Since there is already a 
large number of publications on the topic, here we 
present only a very brief overview of the approach 
of Maximum Likelihood map-making techniques. 

Assuming the simple linear model given by 
Equation [1] the log-likelihood of the data can be 
calculated under the assumption that the noise is 
Gaussian and stationary. The solution is 

logL(d\s) = --{d - AsfN^id - As), (4) 

where N = (n.n*) is the noise covariance matrix in 
the time domain, and . denotes transpose. Maxi- 
mizing the above equation with respect to the map 
parameters s (suppressing the pixel indices here 
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for convenience) leads to the following well known 
estimator: 



a = (A*N- l A)- 1 A t N- l d. 



(5) 



The inverse pixel-pixel covariance matrix of the 
noise in the map is the term in brackets in this 
equation, i.e. 



N- p ) = A t N~ 1 A. 



(6) 



Computation of the solution to Equation [5] is far 
from trivial for most astronomical applications, 
due to the large amount of data, and hence this 
poses a difficult numerical challenge. The noise 
covariance matrix N is a very large matrix of size 
the number of samples squared, which could eas- 
ily be millions, while N pp > may be more reasonable 
in size but has no obvious symmetries, and so is 
still difficult to invert. Nevertheless, we have im- 
plemented a method aimed at finding the Max- 
imum Likelihood solution given by Equation [5] 
when there are a large number of detectors and in 
the presence of possible correlations in the noise 
between different detector timestreams. 

3.2. Implementation 

In the simple case of dealing only with indepen- 
dent noise between detectors, our matrix-inversion 
method is ve ry similar to the MAD CAP method, 
described in Stompor et al. ( 20021 ). However, we 
have developed our new approach to deal effi- 
ciently with multi-detector data in the presence 
of correlated noise between detectors (described 
in detail in Section l3.4p . In this section, we sum- 
marize the basic ideas for the simpler 1-detector 
I-scan case. 

In order to find the Maximum Likelihood so- 
lution of the map (Equation [5]) , we have devel- 
oped two different algorithms. They both allow 
us to solve the linear system N pp is = x, with x = 
A t N~ 1 d. The first approach explicitly computes 
the inverse pixel-pixel covariance matrix N~ p ] , and 
we refer to this as the 'brute- force algorithm'. The 
second approach uses iterations which converge to 
the Maximum Likelihood map without the need 
for computing N~ p ) , and we refer to this as the 
'iterative algorithm'. Both approaches require as 
a first step the computation of the inverse of the 
time-time noise covariance matrix N. 



3.2.1. Inverse noise covariance matrix N tt , 

In practice, even when we have knowledge of 
the statistical properties of the data as described 
by the power spectrum P(u>), the brute- force in- 
version of N is not tractable because of its enor- 
mous size - for a single BLAST detector observing 
for 10 hours at a data-rate of 100 Hz, the matrix 
has approximately 10 13 elements. However, if we 
make the approximation that each data segment is 
"circulant" , meaning that the beginning and the 
end of a segment are connected without disconti- 
nuity and that there are no gaps in the data, then 
the matrix N is also circulant (see Section 13.61 for 
a description of how we treat gaps in the data). 
Circulant matrices are much easier to store and to 
invert. With this approximation the matrix can 
be written 

N tt '=C(\t-t'\), (7) 

where the correlation function C(\t — 1'\) between 
samples t and t' depends only on the separation 
between the two samples. A circulant matrix has 
the property of having a diagonal matrix counter- 
part in Fourier space. 

Let F be the discrete Fourier operator, we have 



N = F^KF, 



(8) 



where denotes transpose conjugation, and A is 
a diagonal matrix whose diagonal is described by 
the power spectrum of the data segment, 



A., 



P(w). 



(9) 



The inverse of the noise covariance matrix is 

TV" 1 = FU^F, (10) 

and because A -1 is a diagonal matrix, A^ _1 is also 
circulant, so that knowledge of only one row is 
enough to describe the entire matrix. Then the 
inverse covariance matrix can be written 



[N -% l= C(\t-t'\), 



with 



C(At) = T~ 



(At), 



(11) 



(12) 



where T~ x represents the inverse Fourier trans- 
form. The inverse of the covariance matrix can 
then be computed directly using the power spec- 
trum of the data. This is a very fast operation 
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(0(n s logn s ), with n s the number of samples), and 
does not require large memory since only one row 
of the matrix is computed. 

The approximation that the data segments are 
ring-shaped or circulant might seem unreasonable, 
but in the end this only has an effect on a small 
fraction of the matrix. For data with large-scale 
correlations (data described by a 1// power spec- 
trum, for instance), the approximation implies an 
assumption that the two edges of the segment are 
very correlated. In reality, there is obviously little 
correlation at long timescales compared to short 
timescales, so extra striping could be introduced 
in the maps if one steps across the two edges of 
the data segment. We have addressed this prob- 
lem in two ways in order to avoid introducing ar- 
tifacts in the final map. In the case where we 
explicitly compute the inverse pixel-pixel covari- 
ance matrix N~Ji (as in the brute-force inversion 
algorithm) , for the estimation of N entering in the 
computation of A~;, we constrain C(At) = for 
At > n s /2, with n s being the number of sam- 
ples in the segment; this is sometimes known as 
"the MADCAP approximation" in the literature. 
It is important to note that this approximation 
cannot be used for the computation of A t N~ 1 d, 
because it is performed partly in Fourier space 
(see Section 13.2. 3|) . Instead we have apodized the 
data d at the edges, and we have removed a low- 
order polynomial (5th order in practice) to reduce 
fluctuations having typical timescales of order (or 
larger than) the data segment (see Section [2.2} . 
This is reasonable, since those scales are not well 
described by a limited number of Fourier modes, 
and hence will always be hard to reconstruct. 

3.2.2. Inverse pixel-pixel covariance matrix N~J; 

The computation of the inverse pixel-pixel co- 
variance matrix, which is described in this section, 
is required only in the brute-force inversion algo- 
rithm, or for an accurate error estimation in the 
pixel domain. 

Since the pointing matrix has only one non- 
zero element per row in our simple model (one 
data sample is associated with a single map pixel), 
the matrix multiplication A t N~ 1 A requires a sin- 
gle loop going across all the non-zero elements of 
A -1 . For most dominant faction of the 

map-making computing time will be devoted to 
this operation. If the data are only correlated 



within a typical length A c , we have the property: 
C(At) ~ for At > A c , and A -1 is a band-diagonal 
matrix (elements separated from the diagonal by 
more than A c are negligible) . The number of ele- 
ments to go through in the loop is of the order of, 
but smaller than 2n s A c , which is hopefully much 
smaller than the size of the matrix itself. 

Unfortunately, if the noise is described by a 
power spectrum of the form (l//) ,a , the correla- 
tion length of the noise is basically of the order 
of the length of the whole data-set. However, the 
amplitude of the correlation is decreasing for very 
long timescales and becomes negligible beyond a 
certain scale. The function C(At) can then be 
artificially set to zero for At > X' c , with X' c cho- 
sen empirically so that the correlation of the noise 
is low enough for scales longer than A c , and also 
that there is very little constraint on the signal 
at those scales. For the specific case of BLAST 
observations, we find that X' c ~ 200 s is a good 
compromise, as illustrated in Figure [TJ 

The noise strongly dominates the signal in 
BLAST observations for scales longer than A c 
(corresponding to frequencies smaller than 5 x 
10~ 3 Hz), since: (1) this frequency is well below 
the knee frequency of the noise power spectrum; 
and (2) there is very little signal at frequencies 
smaller than 5 x 10~ 3 Hz, which is more than ten 
times smaller than the scanning frequency. There- 
fore we have used A^. = 200 s for the computation 
of the noise covariance matrix. The impact of fix- 
ing C(At) — for At > X' c in the initial power 
spectrum is shown in Figure [TJ We can see a tight 
relation between the power spectrum and the in- 
verse covariance matrix, since getting C(At) = 
for At > X' c can be obtained only by modifying 
the power spectrum at frequencies smaller that 
1/A C . 

3.2.3. Computation of A t N~ 1 d 

In the computation of x = A t N~ 1 d, which is 
necessary for both of our algorithms, the multipli- 
cation N~ 1 d is performed in Fourier space where 
the noise covariance matrix is diagonal. We ob- 
tain this vector by dividing the Fourier transform 
of the data by the power spectrum of the noise. 
Another way to represent this is by considering 
that since A -1 is circulant, N~ x d is a convolution 
operation. Assuming that this model holds, the 
resulting data vector d = N~ x d contains whitened 
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Fig. 1. — Left panel: Absolute value of the first row of the inverse covariance matrix given by C{At) (see 
Section f3 . 2 . 1 [) for a typical BLAST observation. The vertical line indicates the value of A^., such that for 
At > X' c , C(At) is set to zero for the computation of N^ 1 . Right panel: Power spectrum of the noise in a 
typical BLAST observation (dotted curve), and after thresholding at frequencies smaller than 5 x 10~ 3 Hz, 
i.e., 1/X' C (solid curve). The dot-dashed curve is obtained by inverting Equation[T2] after forcing the values of 
C(At) to zero for At > X' c , with C obtained from the initial power spectrum of the noise (dotted curve). To 
find the power spectrum which corresponds to the dashed curve the same procedure is applied, but starting 
from the thresholded power spectrum (solid curve) . All the curves are very similar for frequencies larger than 
1/Ap, but the dot-dashed curve begins to diverge at smaller frequencies, while the dashed and solid curves 
lie very close to each other for all frequencies. This shows the tight relation between the power spectrum 
and the inverse covariance matrix. The noticeable peak in the power spectrum is located at the scanning 
frequency of about 4 x 10~ 2 Hz. 



noise. The remaining operation A t d just performs 
the addition of the filtered data sample onto the 
pixels of s (i.e. the map), and hence is very fast. 

In the case of BLAST observations, since we are 
not attempting to recover useful information from 
timescales larger than 200 s, we perform a high- 
pass pre- filtering of the input data d at 5x 10 -3 Hz. 

3.2.4- Matrix inversion algorithm 

In the matrix inversion algorithm, N~J; is di- 
rectly computed, as described in Section 13.2.21 
The next step is to solve the linear system s = 
[N~Ji]~ 1 x. For small maps, in which NZp can be 
stored in memory, we perform a Cholesky decom- 
position. For larger maps, we write the matrix 
to disk and perform an iterative inversion of the 
system using a conjugate gradient method with 
pre-conditioner. 

This algorithm allows one to easily perform 
multiple Monte-Carlo simulations, since the pixel- 
pixel covariance matrix, being independent on the 



data realization, can be computed once (this is as- 
suming that the noise power spectrum is known 
and does not have to be estimated from the data 
themselves), and used for all the realizations of 
simulated data. Another advantage of this ap- 
proach is that it allows for the exact computation 
of the errors in the map, which are given by the 
covariance matrix. 

However, the matrix inversion algorithm is gen- 
erally slower than the iterative algorithm which 
we will discuss next, and can be used only for rel- 
atively small maps; we found that we were limited 
to around 200,000 pixels if the matrix is written 
to disk, or less than 20,000 pixels if the matrix is 
stored in memory. 

3.2.5. Iterative algorithm 

We now present an iterative algorithm based on 
conjugate gradient with pre-conditioner to obtain 
the Maximum Likelihood solution for the map (a 
similar algorithm has been used in lAshdown et al. 
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(120071) 1. Let us rewrite Equation [5l which relates 
the best estimate of the map with the data, after 
multiplying both sides by the pixel-pixel covari- 
ance matrix: 

A t N~ 1 A s = tfN^d. (13) 

If we define Sk as an estimate of the map at iter- 
ation fc, the conjugate gradient method allows to 
solve the linear system by minimizing iteratively 
the following criterion: 

* = r'N-W, (14) 

where 

r^iA'N^Ah-A'N^d). (15) 

This criterion is indeed minimum and equal to zero 
if Sfe is the Maximum Likelihood solution. 

One can interpret \l/ as the weighted variance of 
the difference between two map vectors. The first 
of these vectors, A t N~ 1 A 3k, is the inverse pixel- 
pixel covariance matrix times the current estimate 
of the map, while the second vector, A t N~ 1 d, 
is a map constructed by simply co-adding pre- 
whitened data. We decide that convergence is 
reached when the quantity r*r (which is much eas- 
ier to compute than 'J, and also converges to zero) 
gets smaller than a predefined value. In practice 
the number of iterations required for convergence 
is of the order of 100. 

The conjugate gradient method is not described 
here, since it is a fairly standard numerical tool, 
and the interested reader can find many descrip- 
tions in the literatureQ Instead of describing the 
details, we focus here on aspects which are specific 
to our map-making process. In particular, let us 
describe the computation of A t N~ 1 A St, which is 
the time-consuming part of the optimization and 
has to be performed at each iteration (the compu- 
tation of A t N~ 1 d, also time consuming, but needs 
to be done only once, since none of the parameters 
are changing through the iterations); the other op- 
erations for updating the map at each iteration are 
significantly faster. 

One advantage of this iterative algorithm is 
that the computation of the full pixel-pixel co- 
variance matrix is not required, and the operation 

*e.g. the 1994 article by J.R. Shewchuk: |http : / /www, cs ■ emu ■ edu/ 
~quake-papers/painless-con jugate-gradient .pdf 



A t N~ 1 A Sfe can be done step by step. Indeed, we 
start by computing d = A§k, which is an estimate 
of a "signal" timestream. This operation is equiv- 
alent to scanning over the current estimate of the 
map using the pointing solution. The subsequent 
operation A t N~ 1 d (which should now be famil- 
iar), is carried out in Fourier space, as described 
in Section r3.2.3l (without applying any extra filter- 
ing), and in Section f3T4l for the case of correlated 
noise between detectors. 

This iterative approach is in general much 
faster than the brute-force inversion approach, 
because the most time-consuming operations are 
performed in Fourier space. It also requires less 
memory, since -WZ/ is not explicitly computed. 
Of course, if there are found to be (or known to 
be) non-trivial correlations in N~J; , then it may 
have to be calculated explicitly, hence requiring 
the brute-force approach. However, provided the 
pixel-pixel correlations only involve relatively few 
pixels, it should be possible to calculate a re- 
stricted part of (or perhaps an approximation for) 
N~J in a modified iterative approach. A related 
concept is discussed in Section l3~8l 

3.3. Multi-scan, multi-detector case 

In the previous sub-section, we presented the 
general method for the simple case where only a 
single continuous observation is considered. We 
now describe how we combine observations from 
different detectors at the same wavelength, as well 
as different data segments obtained over differ- 
ent "visits" during the flight, where by "visit" we 
mean a period in the data which starts after a 
sufficiently long gap, or after the observation of a 
different region of the sky. 

For convenience, the data vector d in our model 
(Equation [T]) now contains all the individual data 
segments from different detectors, and also within 
a single channel, concatenated end to end. The 
noise vector n in Equation [T] is defined in a simi- 
lar manner. The matrix A in Equation [1] is then 
the result of stacking individual pointing matrices. 
The Maximum Likelihood solution is also written 
as in Equation^ with N becoming the full covari- 
ance matrix of the noise, including all the channels 
and data chunks. To start with, let us assume that 
there is no correlation of the noise between data 
segments. This is a very good assumption if we 
consider data segments obtained over different vis- 
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its, but is certainly not a good approximation for 
segments obtained simultaneously with different 
channels, since we found a very strong common- 
mode noise between detectors. We will consider 
the simple no-correlation case first, and then in 
the next sub-section fSection l3.4p we will general- 
ize the map-making method to account for noise 
correlations from different detectors. 

In the absence of correlations between data seg- 
ments, the time-time noise covariance matrix N is 
block-diagonal and each block can be inverted sep- 
arately. Defining Ng as the sub-covariance matrix 
for the data in segment number £, and Ag as the 
sub-pointing matrix going from the map to the 
data segment t, the inverse pixel-pixel covariance 
matrix can be written as 

£ 

The computation time for obtaining this matrix is 
proportional to the number of data segments. In 
this simple case the computation of x — A t N^ 1 d 
can be written 

x = J2 A t N r 1 de, (17) 

e 

where the computation of each term is fast and 
can be performed partly in Fourier space. 

3.4. Detector-detector correlated noise 

We now allow for the presence of correlations in 
the noise between different detectors. In the case 
of multiple visits, the noise covariance matrix, N, 
still has null cross-terms for samples from two dif- 
ferent data visits. Therefore, if the data vector is 
sorted by visit, then N is block diagonal and each 
block contains the correlation coefficients between 
all the detectors for the samples within the time 
interval defined as a single visit. Each visit can 
be treated independently, since the sub-matrices 
can be inverted separately, and Equation[in]is still 
valid, but in this case I is the label for the blocks in 
N. In the following, we therefore focus on a single 
visit, and consider observations by all the detec- 
tors; the generalization to multiple visits should 
be clear. 

To simplify the notation, let N denote the noise 
covariance matrix for the visit being considered, 
with d and n being the data and noise vectors, 



respectively, containing the timestream segments 
for all the detectors put end to end, and iVy being 
a block of N of size n s x n s , corresponding to the 
noise correlations between detectors i and j. Let 
us define F the multi-channel Fourier transform 
operator such that 

n = Fn, (18) 

with n containing end-to-end Fourier transforms of 
each data segment. F is a block-diagonal matrix, 
and each block is the Fourier transform operator 
F for one data segment. 

In Fourier space, the noise covariance matrix R 
can be written 

R = FNF^. (19) 

If we consider a single block of the noise covariance 
matrix for detectors i and j, we obtain: 

R ij = FN ij F^. (20) 

Under the assumption that the data are stationary 
and continuous at the edges (see Section l3~2l for a 
discussion), Nij is a circulant matrix, since each 
element [Nij] t t> depends only on the time interval 
\t — t'\. Rij is then a diagonal matrix with the 
diagonal given by the cross-power spectrum of the 
noise between detectors i and j: 

[Rij] UU ' = Pij(w) i£(w = cj') 

= otherwise. (21) 

Here P(ui) is the noise covariance matrix of size 
rid x n<i for a given mode u>, where rid is the to- 
tal number of detectors. The computation of the 
inverse of R is straightforward, since each Fourier 
mode can be treated independently. If P _1 (oj) is 
the inverse noise covariance matrix for mode u>, 
the same relation as in Equation [2T] applies be- 
tween and all P^ 1 . 

From Equation [TH we can calculate the inverse 
covariance matrix of the noise in real space: 

N- 1 = F*R- 1 F. (22) 

Then, a block of ./V -1 between detectors i and j 
can be written 

[N-% =F^[R-%F. (23) 

Because is a diagonal matrix (as discussed 

previously) [iV _1 ]y is circulant, and is related to 
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the inverse of the matrix containing the cross- and 
auto-power spectra of the noise: 

[N- 1 ] ijW =r- 1 {[P- 1 ]ii}V-t). (24) 

From this relation, we can see that in practice TV -1 
is relatively easy to construct, since each of its 
blocks (referring to each pair of detectors) is a cir- 
culant matrix, so only a row of each sub-matrix 
needs to be calculated using the Fast Fourier 
Transform. Finally, in the case the noise covari- 
ance matrix is used for multiplication in real space 
(as in the brute- force algorithm), the same ap- 
proximation described in Section [3.21 is performed 
on each block of N^ 1 , i.e., [-/V -1 ]^' = for 
\t — t'\ > X' c (or for \t - 1'\ > n s /2 if X' c > n s /2). 
The global structure of the final inverse noise co- 
variance matrix is illustrated in Figure O 

In our model of the data (Equation [3]) in which 
all the correlations between detectors are de- 
scribed by a single common-mode, P _1 (tj) can 
be related to the power spectra of the common- 
mode and of the uncorrelated part of the noise 
between detectors: 

P -1 (w) = (a(c(w) t c(w))a* + (n(c i ;) t n(w))) _1 . 

(25) 

This relation can be generalized if multiple 
common-mode components are present: a would 
then be a mixing matrix and (c(w)^c(w)) becomes 
the covariance matrix of these noise components. 

Having computed the inverse noise covariance 
matrix of the timestreams, we can express (us- 
ing Equation [6|) the noise covariance matrix in the 
pixel domain: 

N P P '=Y, A * l N ~% A i> (26) 

y 

where, as before, i and j label the detectors. The 
computation time for N~J; is now proportional to 
the number of detectors squared. The calculation 
of x = A t N^ 1 d is straightforward: 

x = Y,A\[N-\ j d j , (27) 

ij 

and this is fast, since (as already shown) the op- 
eration [N~ 1 ]ijdj is a convolution, which can be 
performed in Fourier space. One can see from 
Equation [24] that x can be expressed directly with 




Fig. 2. — Inverse noise covariance matrix (iV _1 
in the text) for three detectors and for only 10 
minutes of data (corresponding to 60,000 samples 
per detector). The matrix is computed following 
Equation l24i and the approximations described in 
Section [3.4l are applied. The cross- and auto-power 
spectra of the noise, used for the calculation of 
TV -1 , are computed from the data themselves (one 
of the auto-power spectra is shown in Figure [3J 
the cross-power spectra are an unbiased measure 
of the common-mode signal shown in the same 
figure). Each sub-matrix is circulant and corre- 
sponds to a particular pair of detectors. One can 
see that the off-diagonal sub-matrices have ampli- 
tudes of the same order as that of the diagonal 
sub-matrices. This is due to the very high level of 
noise correlation between detectors. 

respect to the cross- and auto-power spectra of the 
noise: 

a r = 2^^ 1 {[P-%( W ).d i («;)}. (28) 

ij 

The formalism presented above can also be gen- 
eralized easily to deal with detectors operating at 
different wavelenths. The map vector s could be 
merging different maps at different wavelengths 
and the noise matrix iV would account for all the 
correlations of the noise between detectors. The 
joint multi-band map-making would be suitable 
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in practice when some contaminations from ther- 
mal fluctuations in the instrument or atmospheric 
emission are present, because they correlate the 
noise at all wavelengths. 

3.5. Noise power spectrum 

The Maximum Likelihood solution for the final 
map depends on the noise power spectra for each 
data-set (through N in Equation 0, which are as- 
sumed to be perfectly known. However, in prac- 
tice the noise power spectra have to be inferred 
from the data themselves and some uncertainties 
are associated with this iterative process. 

In practice a first (approximate) estimate of 
the noise power spectrum can be obtained by re- 
binning the power spectra of each data segment, 
neglecting the contribution of the astrophysical 
signal in the timelines. Indeed, for most of the 
fields observed with BLAST, and in particular 
for blank extragalactic fields (like the ELAIS-N1 
field in BLAST'S flight from Sweden), the noise is 
highly dominant over the sky signal at all frequen- 
cies. However, this is not true for measurements 
of bright regions in the Galactic Plane. There- 
fore, for the first iteration's noise estimate we focus 
on the data taken for about 6 hours while scan- 
ning the deepest extragalactic field (ELAIS-N1) 
and use this to estimate the noise power spectra, 
which then become the noise input for making the 
first set of maps of all of our fields. This approach 
is not entirely satisfactory, since the noise is not 
stationary during the flight - the noise power spec- 
tra are seen to vary over long timescales, although 
they are quite constant within a single visit of each 
field. This non-stationarity appears most obvious 
when there are variations of the scanning strategy 
between different visits (a change of the scanning 
frequency induces a variation of the location of 
peaks in the power spectrum), but can also occur 
due to variations of detector loading or the detec- 
tor bias being changed during the flight. 

Because of the observed non-stationarity of the 
noise, we would like an independent estimate of 
the noise power spectra for each visit for each of 
the observed fields. After starting from a first esti- 
mate of the noise power spectra based on ELAIS- 
Nl data as described above, we adopt an iterative 
approach between maps and noise power spectra. 
At each iteration, the estimated maps are sub- 
tracted from the data, prior to noise power spec- 



trum estimation. We can summarize our proce- 
dure in the following steps: 

• Estimate Pq(lj) from d using a field known 
to have little signal; 

• Compute s from Equation[5](and also Equa- 
tion[M]or[T21 depending on the noise correla- 
tion being considered) using Pq(u>) as input; 

• Estimate P(ui) from d — As; 

• Re-estimate s from Equation [5] using P(u>), 
and iterate on these last two steps until con- 
vergence is achieved. 

We stop iterating when the noise power spectra do 
not vary by more than 1 per thousand from itera- 
tion to iteration (and find that in practice only 3 
to 6 iterations are necessary to reach convergence) . 

We now focus on how we estimate the noise 
power spectra P(u>) in steps 1 and 3. For the 
simplest case, where no correlations are assumed 
between detectors, we simply compute a bin- 
averaged power spectrum for each data segment: 

Pt(q) = — E ( 29 ) 

where, for bin number q, n q — u> max (q)—u) m i a (q) + 
1, £ labels the data segment, and d is the data 
vector from which an estimate of the map has al- 
ready been subtracted. We have chosen logarith- 
mic spacing between bins and an estimate of P{oj) 
for each uj mode is obtained by logarithmic inter- 
polation, which leads to a smooth power spectrum 
estimate. 

In the more complicated case where correla- 
tions between detectors are assumed to be im- 
portant and therefore are not neglected, we must 
estimate, for every iteration, each cross- and auto- 
power spectrum of the data between detectors, 
Pij(u>), which enter into the computation of the 
inverse noise covariance matrix (Equation I24[) . 
Each cross-power spectrum could be directly es- 
timated as in Equation [221 (using d^ and dj u in 
the formulae for detectors i and j), but instead 
we choose to reduce the number of parameters 
to estimate at each step, by assuming that the 
data are described by a common mode between 
detectors plus independent noise (Equation [3]) . In 
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the framework of this model, the expected cross- 
and auto-power spectra depend directly on the 
following parameters: a, the amplitude of the 
common mode in each channel; (c(o>)* • c(w)), the 
power spectrum of the common-mode part; and 
(n(ui)* ■n(ui)) 1 the power spectrum of the noise 
component, which is independent between detec- 
tors. The relation between the model of P(u>) and 
the parameters has been shown in Equation 1251 
These parameters are typically not known a priori, 
and must be measured using the data themselves. 
We use a blind 'component separation' method 
developed for an entirely different prob lem in 
Delabrouille. Cardoso, k. Patanchon (2003). This 
allows us to obtain a single estimate of all the pa- 
rameters described previously, by simultaneously 
using all the observed timestreams of a given field 
for all the detectors in a specified channel (i.e. at 
a single frequency). The method is known to 
be the Maximum Likelihood solution for a Gaus- 
sian and stationary model of both the noise and 
the common-mode. The cross- and auto-power 
spectra P(oj) are then computed following Equa- 
tion [23 using these same estimated parameters. 
Figure [3] shows the estimated noise power spectra 
in a sample of BLAST data for one representative 
detector using three hours of timestreams during 
scans of the ELAIS-N1 field (which is known to 
be essentially devoid of signal). The auto-power 
spectrum is shown, as well as its decomposition in 
terms of the common-mode power spectrum and 
the uncorrelated noise power spectrum. 

Because our estimate of the converged map of 
the sky s is not perfect and contains contribu- 
tions from residual noise, then in subtracting a 
simulated signal timeline from the data to esti- 
mate the noise power spectrum we reintroduce 
some noise to the data which could potentially 
bias our estimate of the noise power spectrum 
(|Ferreira k Jaff dl200d : lHamiltonll2003h . However 
this effect is greatly reduced by the large redun- 
dancy in each pixel of the final maps, as a result of 
the many repeated scans and the large number of 
detectors at each wavelength. The bias can be ne- 
glected to first order for BLAST, since the noise 
level per pixel in the final map is much smaller 
than the noise in individual detector timestreams. 
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Fig. 3. — Total noise power spectrum (solid) 
for one of the 250 /jm detectors, using about 
three hours of data, corresponding to scans of the 
ELAIS-N1 field. The dotted curve correspond to 
the estimated power spectrum of the common- 
mode between detectors at 250 /im, rescaled by 
an amplitude factor for the specific detector being 
considered (a, in the text, where i is the detector 
index), also estimated using the data themselves. 
The dashed curve represents the estimated power 
spectrum of the uncorrelated part of the noise for 
this detector. The common-mode is very strong 
in these data and dominates over the uncorrelated 
noise at frequencies lower than about 0.1 Hz. Most 
of the low frequency noise excess comes from this 
common-mode. The uncorrelated part of the noise 
shows a knee frequency of less than 0.02 Hz, which 
is 5-10 times smaller than for the total noise power 
spectrum. The power of the uncorrelated noise is 
thus reduced by a factor of about 100 at low fre- 
quencies. The excess signal at the scanning fre- 
quency (peak around 0.04 Hz in the power spec- 
trum) is completely common between detectors. 



3.6. Dealing with gaps in the data 

In order to derive the formalism presented so 
far, we have assumed that each data segment is 
stationary and hence consists of a continuous se- 
ries of data points. However, we have seen in Sec- 
tion [2]2] that the BLAST05 data contain multi- 
ple gaps of typical size less than one second. The 
amount of data in these gaps is a few percent of the 
total. In order to reasonably restore the continuity 
of the data we have filled the gaps with random 
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noise, as described in Section [2~2l The data sam- 
ples generated in the gaps are not reprojected into 
the final map but are directed to "dummy" pix- 
els. In principal, the optimal approach would be 
to create one dummy pixel per flagged data sam- 
ple, avoiding the possibility of several simulated 
samples falling on the same pixel (through the re- 
binning of flagged pixels we do not want to intro- 
duce any spurious constraints for the map-making 
process which could arise from adding crossings 
over different time intervals). However, the ap- 
proach of using one dummy pixel per gap sample 
is impractical, because the total number of dummy 
pixels would be excessive for typical BLAST ob- 
servations, and the pixel-pixel covariance matrix 
(which has size the total number of pixels squared) 
would be prohibitive to store and compute. 

We have adopted a simpler approach, which 
does not lead to the mathematically exact solu- 
tion, but comes very close (as has been shown in 
simulations). This consists of rejecting (for the 
computation of NZp) all the elements of N^ t } as- 
sociated with flagged samples. This is equivalent 
to removing from N~, the rows and columns cor- 
responding to the dummy pixels before the inver- 
sion of the matrix, as opposed to after inversion, 
which would be the correct treatment discussed 
above. This approach is also equivalent to as- 
suming null off-diagonal terms for those rows and 
columns. However, such dummy pixels are obvi- 
ously correlated to some degree with the real pix- 
els in the map, and hence this cannot be entirely 
correct. 

Nevertheless, we have verified that this approx- 
imation has a small impact (a few percent only) 
on the final map for most of the angular scales 
which are sampled, although we found some dif- 
ferences at very large angular scales, sizes of the 
order of the map size. For now, we have not put 
much effort into recovering those very large scales 
because they are subject to other effects, as dis- 
cussed in Section 15.2.21 The minimal impact on 
the final map has been verified using pure signal 
simulations and by comparing the results obtained 
between our simple approach and the correct map- 
making solution. This approach works to a high 
degree of accuracy on most scales because the gaps 
are small and do not introduce important discon- 
tinuities in the timestreams. 

We have used this simple approach because it 



gives sufficiently accurate results over the relevant 
angular scales, while being simple and fast to im- 
plement. However, another iterative procedure 
could be adopted, which would lead to the exact 
solution. In this approach, we define two maps. 
The first map "A" is made from only the uncor- 
rupted (i.e. real sky) samples, while the second 
map "B" is obtained from projecting the simu- 
lated (i.e. for gap-filling) samples. The difficulty 
arises in deciding what to do when simulated data 
from different scans or detectors fall in the same 
pixel - one might want the "generated signal" (and 
not necessarily the noise) to be identical in both 
measurements, in order to satisfy the map-making 
hypothesis. If this condition is not satisfied, then 
some artifacts may be introduced into both maps. 
Here is a solution to this problem: 

1. generate a first set of maps "A" and "B" 
after filling the gaps in the timestreams with 
white noise + a linear baseline. 

2. fill the gaps in the data with the best esti- 
mate of the signal in map "A" together with 
white noise + a baseline which is fitted in 
the gap vicinity of the data "minus" signal 
timestream. 

3. recompute the maps "A" and "B" . Step 2 
ensures that the signal is the same for each 
generated sample falling in the same pixel of 
map "B". 

This approach can also be coupled with the pro- 
cedure for estimating the noise power spectra de- 
scribed in Section 13.51 Preliminary results indi- 
cate that this approach works in practice. Detailed 
studies will be presented in a future publication. 

3.7. Pixel constraints 

For some specific observed fields we may have 
strong priors about the sky emission at a given 
location. For instance, we know that over some 
regions the astronomical signal should vary very 
smoothly or should be very weak with respect to 
the noise, at least outside some localized region. 
This is the case in particular when we map bright 
extragalactic sources in order to calibrate the de- 
tectors and estimate the beams; in these cases, 
regions beyond some predefined distance from the 
beam center can be assumed to have null flux (or a 



14 



constant relative flux in the map, since we do not 
have access to the DC level in maps). If we really 
have strong prior knowledge that we are dealing 
with a bright localized region, then we can take a 
further drastic step - we can constrain the map to 
have the same value in some domains of the sky 
by defining a single pixel containing all the data 
samples falling in that region. 

In practice, we define a small box centered at 
the source location and constrain the part of the 
map outside this box to have a constant value. 
This is a very efficient way to remove stripes from 
the map, since the extremities of all the paths 
across the map are re-adjusted. We have used this 
technique to make maps of the isolated calibrators 
observed by BLAST (Truch et al. 2007). 

3.8. Error estimation 

The variance of the noise in each pixel of 
the final map and its correlations are directly 
given by the pixel-pixel covariance matrix N pp i — 
(A t N~ 1 A)~ 1 . This is true given the following 
assumptions: that our model of the data holds, 
in particular that the noise is a purely Gaussian 
random process, which may not be the case in 
practice at low frequencies; and that our estimate 
of the sample-sample noise covariance matrix N tt > 
is accurate enough that the errors do not propa- 
gate significantly into the final map. As already 
mentioned, we never explicitly compute the co- 
variance matrix, but rather its inverse. The direct 
inversion would take a prohibitive computation 
time for most applications. However, to first or- 
der, we can obtain an estimate of the errors by 
inverting the diagonal part of N~ p ) only, neglect- 
ing the off-diagonal terms. This is equivalent to 
assuming that the noise in the final map is white. 
We have checked with the help of simulations that 
this very simple approximation is accurate to bet- 
ter than 10 percent for all our BLAST05 fields, 
even for those with very poor cross-linking. 

Provided that the size of the map is reasonably 
small, so that we are able to explicitly calculate 
Nppi , we can obtain an accurate estimate of the 
errors for a limited (and small) number of pixels in 
the map. The variance for pixel p, and the covari- 
ance with respect to the other pixels of the map, 
can be computed by solving the linear system: 

(np - rip*) = N p < p Up, (30) 



where u p is a unitary vector with a single 1 
for pixel p. If the Cholesky decomposition of 
the matrix has already been performed for the 
map-making procedure, the computation of Equa- 
tion [3D] is relatively fast and hence can be carried 
out for a grid of non-adjacent pixels, for example. 
This can be used to check the validity of the error 
prediction approximation described in the previ- 
ous paragraph. 

3.9. Computational requirements 

For the brute-force inversion algorithm (in 
which the full inverse pixel-pixel covariance ma- 
trix N~, is computed), five minutes of compu- 
tation with a single 3 GHz processor are needed 
to process 2 hours of data from a single detec- 
tor at a rate of 100 Hz. The computational time 
is proportional to the number of samples if this 
is longer than the assumed correlation length of 
the noise in the data (which has been evaluated 
to be A c = 200 s in BLAST05 timestreams) . If 
noise correlations between detectors are also to be 
accounted for, the computational time is propor- 
tional to the square of the number of detectors. 

Most of the computing time is spent on calcu- 
lating N~ p ). Inversion of the linear system to es- 
timate the map s is relatively fast (a few minutes 
to a few hours for maps of several square degrees 
in size). 

For the iterative algorithm, about two minutes 
of computational time is required for two hours of 
data (under the same conditions described above) . 
This assumes that 100 iterations are necessary to 
reach convergence (which is a realistic number for 
most applications), and the algorithm scales with 
?i s logn s . The situation is much better than for 
the brut-force inversion algorithm if correlations 
between detectors are included. In that case, the 
algorithm scales with the square of the number 
of detectors if this exceeds about 40. If there are 
fewer than about 40 detectors then the algorithm 
scales linearly with the number of detectors. As an 
example, if there are 100 detectors, then including 
noise correlations between the detectors increases 
the computational time by a factor of four with 
respect to the "no-correlation" case. The full pro- 
cessing of 10 hours of BLAST05 data at all wave- 
lengths, including detector-detector correlations, 
can be done with a single processor in a few days. 
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In terms of memory, the brute-force inversion 
algorithm requires storage of the full NlS, matrix. 
However, for the iterative algorithm, only vectors 
of the size of the maps need to be kept in memory, 
which is much less demanding. 

4. Simulations for testing SANEPIC 

We now focus on the application of SANEPIC 
to data. Our aim is to develop tests to validate our 
method using simulated BLAST observations. We 
derive conclusions about how well low frequency 
noise in the maps can be reduced, depending on 
observational parameters such as scanning strat- 
egy, and we compare the results obtained with 
those from simpler methods based on filtering the 
data, e.g., common-mode subtraction. In this sec- 
tion, we describe the simulations performed to test 
the SANEPIC method. 

We have generated several different sets of sim- 
ulations of BLAST timestreams. Each set of simu- 
lated timestreams, representing one particular ob- 
served field, is generated for all the BLAST de- 
tectors used for the analysis of real data (132 at 
250 /mi, 78 at 350 /mi and 39 at 500 /mi) and is the 
sum of simulated astrophysical signal, indepen- 
dent noise, and common-mode noise between all 
detectors. The noise is generated randomly with 
Gaussian statistics, given fixed power spectra de- 
rived from real BLAST05 data. Figure [3] shows an 
example of the power spectrum of the noise in the 
data used as input to the simulations for one of 
the fields. The part of the noise which is indepen- 
dent between detectors is generated for every de- 
tector timestream and has a power spectrum well 
described on average by a relatively flat plateau 
for frequencies larger than about 0.05 Hz, and by 
a part proportional to (l//) 2 5 for smaller fre- 
quencies (these characteristics vary slightly from 
detector to detector). Our knowledge of the real 
bolometer noise power spectrum at low frequency 
is limited by the very dominant common mode. 
In simulations, the common- mode noise is gener- 
ated once for all detectors and has a power spec- 
trum very well fit by a power law with an index 
equal to approximately 2.5, together with some 
broad peaks, the largest being at the scanning fre- 
quency (the amplitude of the peak depends on the 
scanning strategy and the observed field). The 
common-mode power spectrum has an amplitude 



such that it reaches the level of the independent 
noise at about 0.3 Hz (see Figured]). The gener- 
ated common-mode timestream is multiplied by 
an amplitude factor which varies from detector to 
detector by ~ 10% and is added to the simulated 
detector timestreams. The amplitude factors used 
for the simulations have been estimated from the 
data themselves. 

In order to represent the astrophysical signal, 
we have simulated simple maps of diffuse emission 
with a power spectrum proportional to fc -3 , as 
for typical Galactic cirrus emission (e.g. Miville- 
Deschenes et al. 2007). Maps are generated fol- 
lowing Gaussian statistics with a resolution of 1", 
much higher than the typical pixel sizes in the fi- 
nal maps, in order to reduce artifacts due to re- 
pixelization. The amplitude of the fluctuations of 
the simulated map is chosen to match the expected 
level of signal in each observed field. The simu- 
lated maps are scanned using BLAST05 pointing, 
and pure signal timestreams are generated for each 
detector. Signal and noise timestreams are added 
at the end of the procedure (but see Section [5] for 
an explanation of why this operation is not always 
carried out). 

We have generated two sets of simulations 
which correspond to two different fields observed 
by BLAST. We selected two fields that were ob- 
served with very different scanning strategies, 
since the performance of the map-making pro- 
cedure is very dependent on scanning strategy; 
this allows us to test SANEPIC in two very differ- 
ent configurations. In the first case the scanning 
was performed mainly in a single direction over a 
short time interval, while in the second case the 
field was observed several times during the flight 
at different scanning angles, to achieve significant 
cross-linking in the map. 

The first data-set uses observations of the Cas 
A supernova remnant emission which comprises 
about 20 minutes of data. BLAST observations of 
this field and derived conclusions will be described 
in detail in Hargrave et al. (in preparation). The 
rectangular region mapped has a size of the order 
of 0.5 deg 2 and was scanned two times back and 
forth over a short time interval. We have gener- 
ated simulations corresponding to all the detectors 
at 250 /mi (we used a total of 132 detectors). 

The second data-set reproduces the obser- 
vations of the intermediate velocity cloud IVC 
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G86.5+59.6 (hereafter 'G86'). Simulations include 
four different visits of the field performed during 
the flight at very different time intervals (rang- 
ing from a few hours to more than a day). Each 
continuous observation segment has a size which 
varies from one to two hours. Two scanning di- 
rections are dominant, which form an angle close 
to 50°. The region covered has a size of about 
2 deg 2 on the sky. Simulations for this field are 
performed specifically for all the 500 /jm detectors 
(41 detectors used). Similar Monte Carlo simula- 
tions at 250 /im would have taken a factor of 10 
longer, while we believe that the conclusion would 
remain unchanged. 

A total of 20 sets of simulations of the observa- 
tions for each field have been performed. For each 
set we vary the realization of the noise and of the 
signal input map. About four hours of computing 
time are needed to create one realization of a full 
set of simulations of G86 with a single processor, 
compared to a few minutes for the Cas A simula- 
tions. This is using the pre-computed full pixel- 
pixel covariance matrix, which was also used to 
analyze the real data. 

5. Results from simulations 

We now present the results obtained with 
SANEPIC applied to the two sets of simulated 
data. In each case, we compare the final map with 
other maps obtained using simpler map-making 
procedures. For these tests we have assumed 
that the noise power spectra are perfectly known, 
rather than estimated separately from each data- 
set; in practice we fix the noise power spectra to 
be the ones from the simulations. We have verified 
that relaxing this constraint has almost negligible 
change on the final maps. 

For these simulated data-sets, we have applied 
some pre-processing of the timestreams before ap- 
plying SANEPIC, just as we do for the real data. 
We have systematically removed a 5th order poly- 
nomial from each timestream segment and weakly 
high-pass filtered the data, as described in Sec- 
tion 12.21 Finally, for the gap-filling we flag the 
simulated data at the same locations as in the real 
data in order to check the influence of the flagging 
procedure in the final maps. 

In the following we have applied SANEPIC in- 
dependently to pure noise timestreams (contain- 



ing independent noise and common-mode noise, 
but without simulated astrophysical signal) and to 
pure signal timestreams. This procedure allows us 
to easily derive conclusions about the noise prop- 
erties in the final maps, as well as about the signal, 
without biasing the results, because SANEPIC is 
a linear method (as shown by Equation [5]) • This 
is only strictly true if the noise power spectrum 
is fixed as done here, and not estimated simul- 
taneously along with the maps. Then, applying 
SANEPIC on pure noise timestreams and pure 
signal timestreams independently and adding the 
two final maps is rigorously equivalent to applying 
SANEPIC on signal plus noise timestreams. This 
has been checked numerically, and we find that the 
difference is consistent with double floating preci- 
sion error. An important consequence of this is 
that the properties of the noise in the final map 
are independent of the signal-to-noise ratio. 

5.1. Case without cross-linking 

5.1.1. Noise-only timestreams 

We first study the maps resulting from the 
noise-only timestreams in the configuration of Cas 
A observations. The chosen pixel size of the map 
is 25" and matches the pixel size of the maps dis- 
cussed in Har grave et al. (in preparation). We 
compare the noise maps obtained from three dif- 
ferent procedures: 

• Case 1: use SANEPIC with the correct 
treatment of the correlated noise. 

• Case 2: use SANEPIC fixing the correlation 
of noise between detectors to zero and fix- 
ing the noise power spectrum for each de- 
tector to the power spectrum of the sum of 
uncorrelated noise and common mode. This 
procedure is very similar to more standard 
map-makers in the literature (e.g., Stompor 
et al. 2002). 

• Case 3: make a simple re-projection of the 
data onto a pixelized map by simply aver- 
aging the data falling in each pixel, after 
having filtered the timestream data with 
the same very weak low-pass filter used 
for SANEPIC. This procedure is sometimes 
called "co-addition" . 
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Fig. 4. — Final maps computed from simulated pure noise timestreams in the configuration of the BLAST05 
Cas A observations, which have a dominant scan direction. From left to right: maps obtained with SANEPIC 
including noise correlations; SANEPIC with no noise correlations included in the model; and simple pixel 
binning (see text for more details). Note the extended dynamic range of the simple co-added map (right 
panel). The maps have a size of about 40' in the cross-scan direction and about 1° along the scan. The pixel 
size is 25". 



Figure |4] shows computed noise maps for one 
of the realizations of the noise in each of the 
three cases. As expected, the map obtained with 
the simple pixel binning approach contains a very 
large amount of low frequency noise, with strong 
striping visible along the scan direction. Residual 
low frequency noise can also be seen in the map 
obtained using SANEPIC without accounting for 
the noise correlations between detectors. We do 
not expect this method to be very efficient, since 
it is very non-optimal in cases (such as this ex- 
ample) where a very large fraction of the noise 
is correlated between detectors. In contrast, the 
noise map obtained with SANEPIC is quite satis- 
factory, showing reduced power at low frequency 
as compared to the previous case. Nevertheless, 
some very weak excess power is seen in the cross- 
scan direction. This is expected, since the map is 
not cross-linked, and very poor constraints can be 
put on the cross-scan directions at low spatial fre- 
quencies (two positions in the map separated by 
more than the size of the array in the cross-scan 
direction are observed far apart in time). 

In order to quantify the level of low frequency 
noise in the maps, we compute the 1-D power spec- 
tra of the maps, averaged over the 20 realizations 
of the simulated data. For the computation of 
power spectra, we take into account only the cen- 



tral part of each map, where the level of redun- 
dancy in the observations is high (we use only the 
highest signal-to- noise region in the maps). To do 
so, we apply an apodized mask to the maps go- 
ing smoothly from at the edges to 1. Figured] 
shows the noise power spectra in the three cases. 
The noise level in the simple re-projection map 
is obviously very poor at all scales. Both of the 
other map-makers reach the white noise level for 
scales smaller than 3' and have excess power at 
larger angular scales. Nevertheless, the gain be- 
tween full SANEPIC and SANEPIC without cor- 
relations is very important at all scales larger than 
about 2' and reaches a maximum value of about 
10 at around 20' angular scales. An interesting 
fact is that the knee frequency of the noise power 
spectrum in the optimal case here corresponds to 
the inverse of the physical scale of the detector 
array in the cross-scan direction (which is of the 
order of 6'). Indeed, there are no observational re- 
dundancies on scales larger than the array in the 
cross-scan direction in the absence of cross-linking 
in the map. Thus the very long timescale 1/ / 
noise present in the timestreams is not efficiently 
removed and re-projects in the final map at large 
angular scales. This effect is also present along the 
scan direction, but with a lower amplitude as the 
map is scanned back and forth. The trend of the 
large angular scale power spectrum of the noise in 



18 





Fig. 5. — One-dimensional power spectra of the 
noise (rebinned in frequency) in the final noise 
maps after map-making in the BLAST05 Cas 
A configuration. Power spectra are averaged over 
the 20 realizations of the simulated data. The 
dashed curve is for the simple re-projection map, 
the dot-dashed curve for SANEPIC with no noise 
correlation between detectors and the triple-dot- 
dashed curve for SANEPIC including a treatment 
of the correlations. The straight line indicates the 
level of white noise in the map predicted by the 
map-making procedure (see Section 13. SJ) , Error 
bars are computed from the dispersion of mea- 
surements among the realizations. For compar- 
ison, the upper dotted curve (decreasing almost 
like a power law at all scales) represents the power 
spectrum of the pure simulated signal in the final 
map. The solid curve represents the power spec- 
trum of the final map obtained with real data us- 
ing SANEPIC, with correlations included. This 
shows the benefit of taking into account correla- 
tions of the noise between detectors in the map- 
making procedure, reducing the noise structure far 
below that of the signal in the map. The real data 
power spectrum shows that the signal dominates 
at all angular scales larger than about 3' and at 
smaller scales we can see that white noise at the 
expected level dominates in the map. The drop of 
power at around a 3' scale is due to the BLAST05 
beam. 

the map just follows the trend of the low frequency 
noise power spectrum in the timestreams. We will 
see in Section \5l2\ that this effect is reduced when 
there are multiple scanning directions in the map. 
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Fig. 6. — Two-dimensional power spectrum of the 
noise maps in the BLAST05 Cas A observational 
configuration obtained with SANEPIC (noise cor- 
relations included) plotted on a logarithmic con- 
trast scale. 

In order to determine the direction in which 
the noise power is strongest in the map, we have 
also computed the 2-dimensional noise power spec- 
trum. The map of the 2-D power spectrum of the 
noise obtained with SANEPIC (noise correlations 
included) is shown in Figure [5] The large bright 
spot around the center corresponds to a relatively 
isotropic component of correlated noise (at least at 
large angular scales). It contains a large fraction of 
the noise power at large angular scales (seen in the 
1-D power spectrum in Figure 0. A smaller, but 
significant fraction of the correlated noise is con- 
centrated in directions perpendicular to the scan 
direction, as can be seen in the figure. As already 
discussed, the reason for this excess power is that 
the noise in the cross-scan direction is poorly con- 
strained. This cross-scan component of the noise 
is significant all the way up to the pixel scale. 

5.1.2. Signal-only timestreams 

We now focus on the signal-only timestream 
simulations. In order to demonstrate the supe- 
rior performance of SANEPIC relative to sim- 
pler methods based on data filtering, we com- 
pare with a map-making method which consists 
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of the following: we first remove the whole ar- 
ray average from each detector timestream and 
then make maps using SANEPIC, assuming no 
correlations between detectors. Removing the ar- 
ray average reduces the signal to almost zero for 
scales larger than the array and so we expect no 
large scale structures to survive in the map. This 
SANEPIC "common-mode subtraction" method 
is still a better procedure than just reprojecting 
the data (after common-mode subtraction) with a 
well chosen filtering to suppress noise drifts (at 
/cut = 0.02 Hz, for instance, since that corre- 
sponds to the knee frequency of the independent 
part of the noise). This latter method is commonly 
used in the submillimeter community and is re- 
ferred to_asJ|sky_remoyal" in reduction of SCUBA 
data (Jenncss ct al. 

Figure [7] shows the input map for one of the sig- 
nal realizations (left panel) , as well as the maps ob- 
tained with SANEPIC (correlations included, cen- 
tral panel) and with the common-mode subtrac- 
tion method (right panel). Results are expected 
to be worse in the second case, because of the ex- 
tra filtering and also because SANEPIC gives less 
weighting to modes at lower frequency which are 
more contaminated by independent noise. We can 
see from Figure [7] that part of the very large scale 
fluctuations with sizes of the order of the map are 
removed using SANEPIC, but apart from those 
very large scales, the input map and the SANEPIC 
map look very similar. More differences can be 
seen in the map obtained with the common-mode 
subtraction method. This is quantified in FigureEJ 
which compares the 1-D power spectra of the two 
output maps, averaged over 20 simulations and 
multiplied by k 3 . Recall that the input spectrum 
varies as k~ 3 and so deviations from a flat line are 
the result of the map-making reconstruction. Note 
that the vertical scale is linear in Figure [H 

With SANEPIC, the power of the reconstructed 
map decreases for scales larger than about 30'. 
There are three reasons for this: the power spec- 
trum is computed over only a small fraction of the 
sky (and for an apodized map), so that structures 
of the order of the map size are never fully re- 
covered; there is weak filtering of the timestreams 
at /cut = 5 x 10~ 3 Hz and through the 5th order 
polynomial subtraction; modes in the maps that 
are very weakly constrained in the map-making 
procedure tend not to be reconstructed through 
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Fig. 8. — 1-D power spectra of signal-only maps 
for the BLAST05 Cas A field reconstructed using 
SANEPIC (solid curve) and using the common- 
mode subtraction method (dot-dashed curve). 
Power spectra are multiplied by fc 3 , displayed on 
a linear scale and averaged over 20 realizations of 
the simulations. The large angular scale behavior 
shows the effectively filtering of each map-making 
procedure, while the drop off at small scales is 
caused by the PSF. 

the matrix inversion procedure, since the matrix 
is very ill-conditioned and numerical problems oc- 
cur. The last two effects are the dominant ones. 
As a result, modes which are preferentially filtered 
are those which lie perpendicular to the scan di- 
rection. 

At angular scales smaller than about 2', the 
power slightly decreases due to the smoothing ef- 
fect of the pixelization. The transfer function at 
those scales is well described by a sine function. 

Turning now to the common-mode subtraction 
method, the power in the map is significantly re- 
duced for scales larger than about 10', and drops 
rapidly to zero. This is because the common-mode 
subtraction removes power on all scales larger than 
the array. On smaller scales, the filtering effect is 
relatively weak, and is reduced when the number 
of detectors increases. 

For these particular simulations the common- 
mode subtraction method (using SANEPIC, but 
with no correlations) does not in fact perform 
very poorly compared to the SANEPIC optimal 
approach. This is because the observed field is 
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Fig. 7. — From left to right: map of simulated signal used as input for one of the realizations of the simulations 
in the configuration of the BLAST05 Cas A field; reconstructed map using SANEPIC accounting for noise 
correlations; reconstructed map using the simple common-mode subtraction method (intensity units here 
are arbitrary). 



small, with a size just a few times bigger than the 
array, and structure at scales smaller than the ar- 
ray size are not strongly affected. This particular 
map is also not cross-linked. However, the situa- 
tion is different for large cross-linked maps like the 
Vulpecula field, as discussed in Section HT2"1 

5.2. Case with cross-linking 

5.2.1. Noise-only timestreams 

We now focus on the set of simulations of the 
G86 field at 500 /im. As in the previous example, 
we first examine the maps resulting from noise- 
only simulated timestreams using three meth- 
ods: optimal SANEPIC (with noise correlations 
taken into account); SANEPIC without consider- 
ing noise correlations between detectors; and the 
simple co-add method. The chosen pixel size for 
the reconstructed maps is 1', which allows for in- 
version of the covariance matrix with a single pro- 
cessor and hence rapid Monte Carlo simulations. 
The conclusions drawn would remain unchanged 
if the pixel size was reduced. 

Figure [9] shows the final noise maps in the three 
cases for one realization of the simulations. The 
simple re-projection map is obviously very stripy 
and would be of little use as an estimate of the 
signal; nevertheless it helps to visualize the direc- 
tions of scanning in the map. We can see two 
main directions covering the central region of the 



map, oriented at about 50° to each other. A third 
scanning direction is also visible but has a much 
smaller weight. The central part of the map is 
the cross-linked region, where we expect the more 
optimal map-making procedures to excel. 

In the map obtained using SANEPIC without 
considering noise correlations between detectors 
(middle panel of Figured]), some large scale noise 
is still visible in the map, even if almost no resid- 
ual striping is apparent in the cross-linked region. 
Indeed, too much weight is given to the large 
timescales in the timestreams (which are basically 
common between all the detectors) as compared 
to the smaller timescales for which there are more 
independent measurements, because the degree of 
correlation of the noise is weaker. The residual 
noise at large angular scales is much weaker in 
the SANEPIC map in which we account for the 
proper correlations of the noise between detectors 
(left panel of Figure [9]). The noise in this map 
looks particularly white. 

The noise level in each simulated map is quan- 
tified in Figure 1101 showing the 1-D power spec- 
trum of residual noise averaged over 20 simula- 
tions, and computed over the cross-linked region 
(which has a diameter of about 100'). While sev- 
eral orders of magnitude are gained in the noise 
power at all scales using the SANEPIC "no cor- 
relation" method as compared to the simple re- 
projection method, accounting for the correlations 
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Fig. 9.— Noise-only simulations (like Figure SJ) for the BLAST05 G86 scanning configuration. The maps 
have a size of about 100' across the short axis and 2.5° across the long axis. The pixel size here is 1'. The 
three panels show the map obtained with SANEPIC including noise correlations (left), SANEPIC with no 
noise correlations in the model (middle), and the simple co-added map (right). Note the different brightness 
scale chosen for the last map due to its larger dynamic range. 



with SANEPIC allows us to further reduce the 
noise power on scales ranging from 20' to the size 
of the map by an additional factor of ~ 5. Toward 
smaller scales, the gain between the SANEPIC 
correlation versus no correlation test cases is still 
very significant down to about 10', where the both 
methods start to approach the white noise level. 
In the optimal map obtained with SANEPIC the 
ratio between the noise power spectrum at large 
scales and the white noise level is around 20, which 
is relatively small. Figure ITTI shows the 2-D power 
spectra of the noise in the SANEPIC map aver- 
aged over 20 realizations. As expected, the large 
scale noise is more important in directions perpen- 
dicular to the scans. 

5.2.2. Signal-only timestreams 

We now focus on signal-only simulations for 
G86. As in the case of the Cas A configura- 
tion, we compare the performance of SANEPIC 
with respect to the simpler common-mode sub- 
traction method. Figure [12] shows the pure signal 
input map for one realization of the simulations 
(left panel) compared with two recovered maps. 
The first output map is obtained with SANEPIC 
including correlations between detectors (middle 
panel), while the second is obtained by subtract- 
ing the common mode between all detectors, fol- 
lowed by applying SANEPIC, but neglecting noise 



correlations between detectors (right panel). We 
see that the very largest scales are not recov- 
ered by SANEPIC. This is because of the weak 
filtering applied to the timestreams and, to a 
greater extent, because of inversion problems with 
SANEPIC on scales of the order of the size of 
the map (these scales are very poorly constrained 
by the map-making procedure). In the common- 
mode subtraction method, only the very largest 
scales are suppressed by the map-making proce- 
dure itself, but the filtering effect is more dramatic 
and extends to much smaller scales. 

The effective filtering is quantified in Fig- 
ure 1131 which shows the power spectra of output 
maps averaged over 20 simulations of pure signal 
timestreams. Again the power spectra are multi- 
plied by fc 3 for comparison purposes. In this case 
SANEPIC works well on scales up to about half a 
degree, above which it fails to recover structure in 
the map; this limit corresponds to scales of about a 
quarter of the map and larger. The filtering effect 
is much more pronounced in the map produced 
with the common-mode subtraction method, be- 
ing strong for all scales above around 14'. This is 
of course expected, since subtracting the average 
of the array strongly reduces the signal on scales 
larger than the array size. Therefore, in order to 
recover large and intermediate scale structure in 
the maps, it is beneficial to use SANEPIC instead 
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Fig. 12. — Signal-only simulations (like Figure [7|) for the BLAST G86 configuration. The three panels 
show the input signal map (left), full SANEP1C (middle) and simple common-mode subtraction followed by 
application of SANEP1C without correlations (right). 



of other methods that are based on simply filtering 
the data. 

5.3. Advantages of cross-linking 

The relative level of residual noise at low spatial 
frequency in the maps is significantly reduced in 
the G86 observational configuration as compared 
to the Cas A configuration. The fundamental dif- 
ference is that the G86 observations contain multi- 
ple (essentially two for most of the data) scanning 
directions, while the Cas A observations are real- 
ized with only two passes across the field in the 
same direction. Multiple scanning directions give 
a huge number of additional constraints for the 
map-making procedure. In particular, large scale 
structures in the map are much better recovered 
in directions parallel to scans, because the noise 
there is smaller. Thus having multiple scanning 
angles allows for recovery of the sky fluctuations 
for all directions, and ends up giving almost no 
weight to the individual loosely constrained cross- 
scan k-modes. 

Differences in the results for maps from these 
two example scanning strategies can be quantified 
in two ways. First of all, for the G86 scanning 
strategy, the transition between white noise and 
"excess" large scale noise in the map occurs at 
a scale around 10', while the same transition oc- 
curs at a scale of around 3' for the Cas A scan- 
ning strategy. Secondly, the ratio between large 
scale noise power and white noise power is larger 



by more than two orders of magnitude for Cas 
A than for G86. On the other hand, some caution 
should be taken to not over-interpret this compar- 
ison, because the pixel size we used is larger for 
the G86 map (1') as compared to the Cas A map 
(25"), and therefore the number of crossings per 
pixel is greater for the G86 map. Nevertheless, this 
simulation exercise has demonstrated that cross- 
linking in the map is extremely beneficial, partic- 
ularly for recovering the large scale structures in 
the map. 

5.4. Map-making transfer function 

When carrying out a complex data processing 
procedure, it is important to check whether the 
results are biased in any way. We have found that 
the transfer function of the map-making proce- 
dure, defined as the ratio between the amplitude 
of fluctuations in the output pure signal map rela- 
tive to the input map, is not always exactly unity, 
even for intermediate and small angular scales in 
the map. For example, in the Cas A configuration 
at 250 /im the fluctuation amplitudes in the final 
map are reduced by 3% on average as compared 
to the input map, almost uniformly across spatial 
frequencies and directions. This global discrep- 
ancy reaches the level of 9% at 500 /im. Moreover, 
it is also present for the G86 configuration. We 
believe that this reduction is due to the fact that 
the pixel-pixel covariance is ill-conditioned and nu- 
merical imprecision occurs in the matrix inversion. 
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Fig. 10. — Output power spectrum comparison of 
noise-only simulations after map-making (like Fig- 
ure [5]), for the BLAST05 G86 scanning configura- 
tion. These power spectra are computed only in 
the cross-linked region of the maps, which forms a 
large disk of about 100' diameter and can be easily 
identified in Figure [31 The dashed curve is for the 
simple re-projection map (right panel of Figure [9]), 
the dot-dashed curve is for SANEPIC without con- 
sideration of noise correlations between detectors 
(middle panel of Figure [9]), and the triple-dot- 
dashed curve is for SANEPIC including the corre- 
lation treatment (left panel of Figure [9]). The hor- 
izontal line indicates the level of white noise in the 
map predicted by the map-making procedure. Er- 
ror bars are estimated from the dispersion among 
measurements for all the realizations. Residual 
low frequency noise in the optimal SANEPIC map 
is very low, thanks to the multiple scanning direc- 
tions of this field. The situation is much better 
than for the Cas A observational configuration, 
which had a single scan direction (Figure E|) . 

We find that the bias tends to be smaller when 
the number of detectors is larger and also when 
the number of constraints increases, like when 
we have multiple scanning directions, or when we 
map isolated bright sources (presented in Truch et 
al. 2007) and constrain the data outside a defined 
region to have a constant flux (see Section 13.71 for 
details of this procedure) . Since this bias can be 
estimated using simulations, it is straightforward 
to correct for. We have found that it is not always 
important, e.g. for the large Galactic map in the 
Vulpecula region analyzed in Chapin et al. (2007) 
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Fig. 11. — 2-D power spectrum of a noise-only 
simulation reduced using SANEPIC (like Figure^ 
for the BLAST05 G86 scanning configuration. 
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Fig. 13. — Power spectrum comparison of the 
SANEPIC map (solid) versus the simple common- 
mode removal map (dot-dashed) for signal-only 
simulated data, as in Figure but for the 
BLAST05 G86 scanning configuration. Power 
spectra are multiplied by fc 3 , so that a flat line 
would indicate no filtering. The drop off at small 
angular scales is due to the BLAST05 PSF. 

the bias is negligible. 
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Fig. 14. — Map of the Cas A supernova rem- 
nant at 250 (xm made from BLAST05 data using 
SANEPIC including noise correlations between 
detectors. The map is represented in Galactic co- 
ordinates with 25" pixels and has a size of about 



0.5deg z 



6. Application to real BLAST05 data 

Now that we have looked at signal-only and 
noise-only simulations, we now turn to real 
BLAST data. In this section, we show maps of 
two example fields from the BLAST05 data which 
have been obtained using SANEPIC. 

6.1. Cassiopeia A 

Figure [14] shows the map obtained from the 
observations of the Cas A field at 250 /im using 
SANEPIC including full consideration of the noise 
correlated between detectors. Detailed analysis of 
the maps at the three wavelengths is described in 
Hargrave et al. (in preparation). The properties 
of the noise and the transfer function of the sig- 
nal in the map have been studied in detail from 
Monte Carlo simulations in Section 15.11 Results 
from such simulations are used to characterize the 
map. The power spectrum of this map has been 
compared to results from simulations in Figure [5] 
The map structures are relatively smooth, due to 
the BLAST05 point spread function, which has 
a width of the order of 3', causing the drop in 



the 1-D power spectrum below those spatial scales. 
The signal clearly dominates over noise on angular 
scales larger than about 3', and the diffuse struc- 
ture should be reliable up to a large fraction of the 
overall map size. 

6.2. Vulpecula region 

Another field observed during the BLAST05 
campaign is centered in the Galactic Plane close 
to the open cluster NGC 6823 in the constella- 
tion of Vulpecula. The region mapped has a size 
of about 4 deg 2 and was chosen for its high-mass 
star formation activity. Compl ete analysis of this 
observed field is presented in lChapin et al.l (|2OO70 . 

A few hours of these data were taken at differ- 
ent time intervals during the flight. By design, 
this field has been observed with very different 
scanning directions, and is therefore it should be 
possible to recover diffuse large scale structures. 

The map of the observed region at 250 fim 
obtained with SANEPIC is shown in Figure [15] 
(left panel). For comparison (right panel), we 
have computed another map using a much sim- 
pler method which consists of removing the ar- 
ray average signal at each timestep for all of the 
timestreams and reprojecting the data onto the 
map after filtering. This is like the "sky removal" 
procedure often carried out for ground-based sub- 
millimeter data (see also Section |5.1.2[) . One can 
see that it suppresses almost all the diffuse struc- 
ture in the map. 

No residual striping is visible in the map ob- 
tained with SANEPIC (left panel of Figure HSJ), 
mainly due to the presence of multiple scanning di- 
rections for this field. Large scale structures in the 
map are successfully recovered with SANEPIC, as 
can easily be seen by comparing with the right 
panel of Figure [T5l This recovery applies to scales 
which are significantly larger than the array size. 
However, the resulting effective filtering after ap- 
plying the "array average subtraction" method 
induces negative signals near bright sources in 
the map, while no such filtering effect is seen 
in the SANEPIC map (except perhaps near the 
edges of the map). This shows that optimal map- 
making methods (in the sense of least squares) 
like SANEPIC are better suited to recover point 
sources in the maps as well as diffuse structures. 
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Fig. 15. — Maps of a Galactic Plane region near NGC 6823 in the Vulpecula constellation derived from 
BLAST05 observations at 250 ^m. The two maps are obtained using two different methods: SANEP1C 
(left panel); and simple reprojection after removing the array average signal at each timestcp from all the 
timcstreams, together with filtering (right panel). Maps arc presented in equatorial coordinates with 15" 
pixels. The region mapped has a size of about 4deg 2 . Numerous point sources have been identified in the 
field (see Chapin et al. 2007), their tell-tale shape in the map resulting from the PSF of the BLAST05 optics 
(see Truch et al. 2007). 



7. Conclusions 

Large format detector arrays operating at far- 
IR and submillimetcr wavelengths are becoming 
the norm, rather than the exception. Ground- 
based instruments are plagued with common- 
mode emission arising from the Earth's atmo- 
sphere. And as we have found with BLAST, the 
same applies to high frequency balloon-borne in- 
struments, where we see correlated noise from 
thermal as well as atmospheric effects. There 
is an expectation that even upcoming satellites 
might be faced with similar issues, because of ther- 
mal variations in the spacecraft, for example. In 
general, we expect that correlated noise between 
detectors will be a major issue which all such ex- 
periments have to deal with, and we expect that 
the SANEPIC approach, which we have described 
here, will be widely applicable. Indeed, there 
is evidence from existing arrays (e.g., SHARC-II 
and AzTEC) that once there are many detectors, 
there are multiple correlations between sub-sets of 



the detectors, as well as an overall common-mode 
term. Consequently, one sees correlations between 
contiguous blocks of detectors on the array, or sets 
of detectors which share amplifiers or are other- 
wise coupled through the electronics. Provided 
that these correlations can be investigated and 
their behavior modelled, it is straightforward to 
extend the SANEPIC approach to deal with sev- 
eral distinct sources of correlated noise. Hence 
we expect the SANEPIC approach to be appli- 
cable to future instruments such as SCUBA-2, 
SPIRE, ACT, Planck HFI and others. There are 
also many experiments being planned which use 
large detector arrays to perform sensitive polariza- 
tion measurements, and we see no reason why the 
SANEPIC approach could not also be extended 
to polarimetry. 
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